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We solve numerically the Einstein-Klein-Gordon system with spherical symmetry, for a massive 
real scalar field endowed with a quartic self-interaction potential, and obtain the so-called $■*- 
oscillatons which is the short name for oscillating soliton stars. We analyze numerically the stability 
of such oscillatons, and study the influence of the quartic potential on the behavior of both, the 
stable (S-oscillatons) and unstable (U-oscillatons) cases under small and strong radial perturbations. 
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I. INTRODUCTION 

Oscillatons are non-singular and asymptotically flat 
solutions of the Einstcin-Klein-Gordon (EKG) equations, 
in which both metric and scalar field are fully time- 
dependent [ll-Q- They can be considered as gravita- 
tionally bounded objects made of a real scalar field in 
the classical regime, and were first discovered by Seidel 
and Suen[l|. Oscillatons should be distinguished from 
their complex counterparts, the so-called boson stars, for 
which the space-time geometry remains static [GHSJ. In 
the literature we can find other works descri bing these 
bound objects and their stability properties |lC)l - [lq | . Also 
there have been many works about numerical evolution 
of scalar fiels in different cosmological and astrophysical 
contextJlMi]. 

Ref. presented an exhaustive study of oscillatons for 
the simplest case of a massive, non-self interacting, scalar 
field. It is shown in there that, as in the case of boson 
stars, oscillatons can be classified into stable (S-branch) 
and unstable (U-branch) configurations. S-oscillatons are 
stable configurations under small radial perturbations, 
and they typically migrate to other S-profiles if strongly 
perturbed. On the other hand, U-oscillatons are intrin- 
sically unstable: they migrate to the S-branch if their 
mass is moderate, but they may collapse into black holes 
if their mass is large enough. 

The case of oscillatons with a self-interaction term in 
the scalar field potential remains without study. This is 
in contrast to the case of boson stars, whose properties 
are well known even in the self- interacting case|8l[9l.[20|. 
It is the main purpose of this paper to cover this omis- 
sion, and to study the properties of self-interacting oscil- 
latons. For that we make use of techniques and methods 
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that parallel those used in boson star studies and for 
non- interacting oscillatons |3l [20l. [2l| . which will ease the 
comparison between the two types of scalar objects. 

A summary of the paper is as follows. In Sec. |lll we 
present the equations motion of the self-interacing scalar 
field. In Sec. IIIII we show the resuls of the equilibrium 
configuration and their representative properties, for dif- 
ferent values of the self-interaction parameter. In Sec. lIVI 
we present the results of the dynamical evolution of oscil- 
latons, and the separate analysis of (stable) S-oscillatons 
and (unstable) U-oscillatons. Finally, Sec. |V] is devoted 
to overall conclusions. 



II. MATHEMATICAL BACKGROUND 



The action that describes our self-gravitating system 
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where V{^) is the scalar potential, 

y(<l>) = im|$2 + iA^, (2) 

771$ denotes the mass of the scalar field, and A is the 
quartic interaction parameter. We shall be interested in 
the spherically symmetric case for which the metric is 
written as 

ds^ = -a^{t.x)dt'^ + a^{t,x)dr^ +r'^{de'^ +sm'^ediP'^), 

(3) 
where a{t, r) is the radial function, and a{t, r) is the lapse 
function. 



A. Evolution Equations 

The energy-momentum tensor for the scalar field 
^{t, r) endowed with a scalar field potential V{^) is de- 



fined as 
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We can identify tlie Klein-Gordon (KG) equation witli 
the conservation equations of the scalar field energy- 
momentum tensor, 



T^"".^ = $" 
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where D = {^/ \/—g)dti{yJ—gd^) is the d'Alembertian op- 
erator. Following Ref. [3|, we introduce the first order 
variables 



* = $' , n = a^/a , 



(6) 



in order to write appropriate evolution equations for the 
scalar fields. Hereafter, a prime denotes derivative with 
respect to r, and a dot denotes derivative with respect to 
t. The KG equation is then represented by the following 
set of first order differential equations 
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As for the metric functions, we will use the so-called 
Hamiltonian constraint for the radial function, 
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and the polar-areal slicing condition for the lapse func- 
tion, 
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For numerical purposes, in all the above equations, we 
have introduced the dimensionless quantities r = a;/m$, 
t -J> t/m,^, $ -J> <&/y^i and A = A/(87rG'm|). Eqs. (O 
|9]) are the total set of evolution equations we shall use to 
explore the properties of self-interacting oscillatons. 



III. EQUILIBRIUM CONFIGURATION 



Before we evolve the EKG equations, we 
attention to equilibrium configurations of 
interacting oscillatons. For that, we consider 
tarily a different set of variables: A(t, x) = 
C(t,x) = a} jo?. In terms of these variables, 
stein equations, Gq^ — 8ttT^^, read 

xA^^' , 
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whereas the KG is the single equation 

x 
-A($-|-A$3). 
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The above equations will be solved by introducing in 
Eqs. ([TU]) and (|lll) . the following Fourier expansions for 
the scalar field and the metric functions. 



$(i,x) = ^ 02j_i(a;)cos((2j- l)wt) (12a) 
A(i,x) = ^ Aj{x)cos{2JLot) (12b) 

C(i,x) = ^ Cj(x)cos(2ja;i) (12c) 

where uj is called the fundamental frequency, and Jmax is 
the mode at which the Fourier series are truncated. We 
define C = CJl"^, with D, = Lo/ra,^. 

After some lengthy but straightforward algebra, we ob- 
tain first-order differential equations for the coefficients 
Ao, C2j, $2j-i and *2j-i, where *2i-i = *2i-i- This 
system of equations can be solve with the aid of an in- 
tegrator, like odeint, that uses an adaptive step-size con- 
trol which allows us to achieve some predetermined ac- 
curacy in the solution with minimum computational ef- 
fort (see ^2|)- We also use a second order Runge-Kutta 
method (RK2) in order to verify the convergence of our 
code, since this method uses a fix step-size. Once the 
profiles of Aq and ^2j-i are known, we will get the sub- 
sequent values of ^2j, for j ^ 0, from the algebraic re- 
striction (llOap . 



A. Boundary Conditions 

Appropriate boundary conditions can be obtained by 
analyzing the regularity of the functions A{t, x), C{y, x), 
and $(i, x) at x — ?► 0. First, we realize that A{t, x = 0) = 
1, which implies Ao(a; = 0) = 1, and Aj{x = 0) = for 
j > 1. Second, we find that $'(i,x = 0) = 0, so that 
</,;.(a: = 0)=OVj. 

From asymptotic flatness, we know that A{t, x — > 
oo) = 1, which leads to Af){x — > cx)) = 1, and Aj{x — >■ 
oo) = for J > 1. The scalar field must vanish when 
X — >■ oo, and this implies 4>j(x — >■ oo) = Vj. Likewise, 
we expect C{t,x — > oo) = fl^, and then Co(x — > oo) = 
il^, whereas Cj(a; — > oo) = for j > I. Actually, the fun- 
damental frequency fl is always an output value obtained 
after the solution of Eqs. (jlObl) . ([TUc)) . and fTTj) . 

Note that there are not boundary conditions for Cj{x) 
and 4>j{x) at X = 0. We need to find these values in 
such manner that all boundary conditions at a; ^^ oo 
are satisfied. For each solution, we will fix the value of 
the first scalar Fourier coefficient (j)i{0), and adjust the 



rest of the coefficients <^i<2 and Cj>q until we satisfy 
the boundary conditions [2]; we use a non-Hnear shooting 
method to find these values. Thus, we obtain a set of 
eigenvalues for each value (t>i{0), which in turn defines 
uniquely a given configuration. 



Numerical Results 



We truncate the Fourier expansions ()12|) at Jmax = 2. 
Then, according to what we said about boundary condi- 
tions, we need to calculate the value of the coefficients 
Co, C2, C4, and 03, for each value of the A parameter. 
Table H] contains the values obtained once the boundary 



conditions were satisfied at a;„ 
01(0) =0.28. 



30 for the fixed value 





A =0.0 


A =1.0 


A =2.0 


A =3.0 


Co(0) 

C2(0) 
C4(0) 
-^3(0) 


8.98e-01 
6.90e-05 
2.66e-06 
5.77e-07 


8.91e-01 
1.16e-05 
7.32e-06 
3.69e-06 


8.88e-01 
1.36e-05 
9.11e-06 
6.85e-06 


8.84e-01 
3.54e-05 
1.68e-05 
6.57e-07 



TABLE I. Resulting values of the Fourier coefficients Co, 
C2, C4, and (jf>3, see Eqs. ([12]), after the integration of the 
equations of motion (|10|) : in all cases, <j>i{0) = 0.28, and 
Xmax = 30. For completeness, the first column shows the 
usual, non-interacting, massive case; it can be verified that 
the values obtained coincide with those reported in|3|. 

In contrast to the non-interacting case, the increasing 
of A complicates the solution of the equations of motion 
as we consider larger values of 01. As a consequence, 
in Fig. [3] we observe that, for A = 2, 3, we were able to 
compute the values of masses and frequencies for scalar 
field values up to 0i(O) = 0.6, whereas for A = 0, 1 we 
could compute up to 0i = 0.7. In the same way, as the 
values of A and of 0i (0) are increased, we must decrease 
the size of the numerical domain, represented by Xmax- 

We present in Fig. [1] the results for the metric func- 
tions grr, gtt, and the scalar field $, for configurations 
with 01 (0) = 0.2 and t = 0. This functions were calcu- 
lated for different values of the parameter A. Note that 
the metric functions are non-singular at the origin and 
asymptotically flat for large x. 

Following [2], the total mass Mt is calculated by the 
so-called Schwarzschild mass, 
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and the fundamental frequency f2 is obtained from 



(13) 
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In Fig. [5J we show the values obtained for the total 
masses and fundamental frequencies as functions of the 
central value of the flrst scalar field (Fourier) coefficient 
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FIG. 1. Metric functions (Top) Qr,-, (Middle) gtt, and (Bot- 
tom) the scalar field $ of equilibrium configurations corre- 
sponding to different values of the quartic parameter A. For 
all examples, we set = 0.2, Xmax = 30. 



(14) 01 (0) and for different values of A. We can see that the 
mass of a given configuration increases for larger values 
of A increases, whereas the fundamental frequency J7$ 
decreases. 

In Table ini we show the values of the fundamental fre- 



quency $7$ and of the critical masses ilf$c for each value 
of A. The critical mass is defined as the maximum mass 
attained by an equilibrium configuration; this maximum 
can be seen in Fig. [5] for each value of A. 



A 


Aig>c{mpi/mg>) 


<^lc(0) 


Q^ 


0.0 


0.605 


0.48 


0.862 


1.0 


0.694 


0.47 


0.851 


2.0 


0.770 


0.44 


0.850 


3.0 


0.854 


0.39 


0.848 



TABLE II. Resulting values of the critical mass AI,s,c for differ- 
ent A's. As expected, the mass of the equilibrium configura- 
tions increase for larger values of A, whereas the fundamental 
frequency decreases. There is also a shift to lower values of 
the critical value 0ic(O) as A increases. 
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FIG. 2. (Top) M<i and (Bottom) Q^ as functions of the central 
value of the Fourier coefficient </>i(0) for different values of 
A. We can observe in that there is a shift to the left of the 
critical value (j}ic{0), which corresponds to that equilibrium 
configuration with a critical mass M^c, as A increases. 

For completeness, we show in Fig. [3l for different equi- 
librium configurations, the critical mass M^c as a func- 
tion of the radius Rmax at which the metric radial func- 
tion reaches its maximum value at i = 0. 



FIG. 3. Typical graph of Mq, in terms of the maximum ra- 
dius Rmax , the radius at which the radial metric function Qrr 
reaches its maximum value, see Fig. [l] for different values of 
A. The configurations become more massive and more com- 
pact for larger values of the self-interacting parameter. 



IV. NUMERICAL EVOLUTION OF 
SELF-INTERACTING OSCILLATONS 

We proceed now to the evolution of the EKG equations 
using the method of lines. For the time integrations, 
we use a fourth order Runge-Kutta method (RK4), and 
second order centered differences to discreticize spatial 
derivatives. A second-order Runge-Kutta method (RK2) 
is used for the spatial integration of the metric functions 
at each time level. Special care is needed for Eq. (|7b|) . 
otherwise its discretization would not be second-order ac- 
curate because of the presence of the factor l/x^ in the 
principal part. The transformation that guarantees the 
accuracy of our simulations at the origin is 



n 



d fx^a"^ 



dx^ 



aa($ + A$^). 



(15) 



A. Boundary Conditions 



To properly account for the origin at cc = 0, we use 
the fictitious point xq — — Aa;/2 as in l3|, and take a 
spatial grid of the form x^ = (i — l/2)Aa;. From the 
evolution equation of $, Eq. ([Tal) . we can see that it not 
necessary to apply boundary condition for $, because the 
evolution equation can be integrated all the way from the 
boundary point xo = —A/2 up to the outer boundary 
point. This is possible since the evolution equation ([7a|) 
does not have spatial derivatives. Likewise, we use the 
fictitious point to impose appropriate parity conditions 
on the scalar field variables: 11 is even, and * is odd. 

At the outer numerical boundary, we assume that 11 
behaves as an outgoing wave pulse of the form 



n = u{x — t)/ X , 



(16) 



where u is an arbitrary function. 
Eq. p^ becomes 



In differential form, 



djl + ftn + U/x = . 



(17) 



By using finite difference in Eq. (fT7|) . we can solve it to 
find the unknown boundary value at the new time level. 
Because 11 behaves as an outgoing wave at the boundary, 
so does $. Then, the outgoing wave boundary condition 
applied to $ implies that at the outer boundary 



* 



-n-^/x. 



(18) 



We used this expression to obtain boundary values for ^ 
after the calculation of those of $ and 11. 

Boundary conditions for the metric functions are the 
following. Local flatness at the origin implies that a{x — 
0) = 1 and dxa{x = 0) = 0, and these two conditions to- 
gether imply that a(a;o) = a[xi) = 1 + 0{/S.xY . We use 
these two values at the first and second grid points to in- 
tegrate the second order Hamiltonian constrain outwards 
to obtain a{x). 

As for the lapse function, we impose a — 1/a as an 
outer boundary condition; this is because in vacuum our 
slicing condition implies that we are in Schwarzschild co- 
ordinates. Then, we are assuming that our boundary 
conditions are sufficiently far away as to be always in 
vacuum. Finally, the slicing condition is integrated in- 
wards to obtain the full profile of a{x). 



B. Initial Conditions 

For most of our numerical experiments, we consider 
as initial conditions the equilibrium configurations cal- 
culated in Sec. IIIII hence, 

$(t = 0,x)= ^ 0,(a;), (19) 

./max 

^(t = 0,x) = $'(0,a;)= ^ ^^-(x), (20) 

n(i = 0,x) = $(d,a;) = 0. (21) 

Our interest resides mainly in the equilibrium properties 
of these configurations, and it is first necessary to develop 
different analysis techniques for their study. In this we 
will include, to begin with, initial configurations that are 
a deformation of the original equilibrium ones. There 
is a practical reason behind: equilibrium configurations 
are easy to evolve, whereas general configurations may 
be very difficult to follow as they develop in time. 



A M^c{mpi/m,j,) 


<^lc(0) 


0.0 0.599 


0.47 


1.0 0.686 


0.46 


2.0 0.767 


0.42 


3.0 0.843 


0.39 



TABLE III. Resulting values of the critical mass Af$c for the 
numerical evolution of equilibrium oscillatons using different 
A's, see also Table HTl and text below for details. 



For instance, the critical mass M$c of self-interacting os- 
cillatons, for different values of A, are shown in Table Hill 
The total mass was obtained through the integration of 
the energy density calculated from the time-time compo- 
nent of the scalar field stress-energy tensor, p$ = — T^o- 
Because of the spherical symmetry, this method gives 
identical results to the Schwarzschild mass in Eq. (|13p . 
We notice a little discrepancy with respect to the values 
reported in Sec. IIIII this is because of the evolution code 
is more accurate and does not depend upon the trun- 
cation we had to impose upon the Fourier expansions 
indni. 



1. Code tests 

To revise the accuracy of our numerical calculations, 
we monitored the momentum constraint (llOap . 



P := a^t a;a*n = 0, 



(22) 



and calculated the L2-noria. of the value of /3 across the 
grid as a function of time using three different spatial 
grid sizes, see Fig. 01 It can be seen that the accuracy 
of the numerical code is not altered by the values of the 
quartic parameter A. 
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C. Numerical Results 

We start by calculating some properties of the equi- 
librium configurations found now by the evolution code. 



FIG. 4. I/2-norm of the momentum constraint (I10a|l for three 
different resolutions dx — 0.01 (red), dx = 0.02 (green), 0.04 
(blue) and A = 3. It can be seen that the accuracy of the 
code is kept under control along the runs. 



On the other hand, Fig. [5] shows the mass-differential 
function x^ p for a run with a numerical boundary located 
at Xmax = 100, and for A = 3. We can see that the 
plots are practically the same up to x = 20, but differ 
one from each other at the outer parts of the numerical 
domain. These discrepancies, of the order oi x^ p < 10^'^, 
arise because part of scalar field has been reflected from 
the outer boundary. The error is not significant and we 
were able to keep it under control during the numerical 
evolution. 



■x 10 




FIG. 5. The function x^p in the case A = 3. The small 
discrepancies seen at the outer parts of the numerical domain, 
of the order of x^p < 10~^, appear because of the boundary 
condition used, which is, in the strict sense, only appropriate 
for the massless case. However, the numerical error is kept 
under control all along the runs. 



Fig. ini shows the evolution of the total integrated mass 
for each time level. We can observe a small adjustment 
of the original mass, because there is a small ejection of 
scalar field at the beginning of the run. We can notice a 
steady decay of the mass, which is evidence of an intrinsic 
dissipation of our numeric code; it can be reduced by tak- 
ing a finer spatial grid, and then it cannot be considered 
an intrinsic decay of the oscillatons. 



2. S-branch and Quasi-normal Modes 

We present here numerical evidence that the charac- 
terization of stable, S-branch, and unstable, U-branch 
oscillatons is preserved in the case of the inclusion of a 
quartic interaction, much in the same manner as in the 
case of non self-interacting oscillatons and boson stars. 

To begin with, we show in Fig. [7] the evolution of 
the maximum value of the radial metric function g^x ~ 
a^(i, x), corresponding to an equilibrium configuration 
with (j>i{Q) — 0.28. The outer boundary of the numer- 
ical domain was set at Xmax = 100, and the run was 
followed up to a time t — 5000. It can be noticed that 
the configuration maintains the same oscillatory pattern 




5000 



FIG. 6. Typical case of mass dissipation in the numerical 
evolution of equilibrium configurations. It is caused by the 
numerical error in the evolution, but we confirmed in our 
numerical experiments that can be reduced by considering a 
finer spatial grid. 



at all times and for all values of A. This is an evidence 
of the stability of the oscillations in response to small ra- 
dial perturbations. In this case, the perturbations come 
from the truncation of the Fourier series (IT^ and the 
discretization error of the numerical solutions. 

Fig. [8] presents the Fourier transform of the oscilla- 
tions shown in Fig. [TJ we observe that the maximum of 
the grr oscillates periodically with two distinctive time 
scales. The short-period oscillation (large frequency) 
corresponds to the fundamental frequency used in the 
Fourier expansions (|12p . The large-period oscillation 
(small frequency) is an overall vibration of the configu- 
ration that we identify as the characteristic quasi-normal 
modes of the oscillatons. 

Following the analysis done in [J], we calculate the 
power spectrum of the evolution for the entire S-branch; 
we show in Fig. IHlthe quasi-normal frequency / as a func- 
tion of the total mass Mt- This kind of plots has proved 
useful in the analysis of the evolution of general scalar 
field configurations, see|3l.[20l.[2]|. 



On the other hand, Fig.[TU]presents the evolution of the 
total mass Mt and Rmax {Umax is the value where the 
radial metric function g^ reaches its maximum value) as 
compared with the values corresponding to equilibrium 
configurations. All the cases correspond to the S-branch. 

It is observed that slightly perturbed S-oscillatons are 
not migrating to another S-oscillaton profile, but rather 
they oscillate with a small amplitude around the original 
equilibrium configuration. It is then concluded that S- 
oscillatons are stable under small radial perturbations, 
and that the frequencies show in Fig. [5] are their intrinsic 
quasi-normal modes. 
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FIG. 7. Maximum value of the radial metric function 
grr{t,x) — a^{t,x), corresponding for the initial configuration 
4>i{0) = 0.28 with Ax = 0.01 and At = 0.005, at x^^ax = 100 
and the evolution is shown to t = 5000m~^ for A = (Top) 
and 3 (Bottom). 



3. Perturbed S-oscillatons 

We now turn our attention to the evolution of some 
strongly-perturbed S-oscillatons. We use a Gaussian pro- 
file as a perturbation applied to the original equilibrium 
configurations; see Fig. [11] for an example of this pertur- 
bation in the case of an S-oscillaton with 0i(O) = 0.1, 
so that its mass is increased by 40%. The purpose is 
to analyze whether S-oscillatons are stable under strong 
perturbations, and the conditions to be met for the col- 
lapse into a black hole. 

First, we increased the original mass of the equilibrium 
configuration (see Fig. [5]) with 0i(O) = 0.1 by 40% and 
60% for all A values. The initial masses for A = 1 arc 
Mi = 0.5748 (40%) and Nh = 0.6816 (60%), for A = 1 
are M, = 0.6075 (40%) and Mi = 0.7203 (60%), for A = 2 
are M, = 0.6423 (40%) and Mi = 0.7615 (60%), and for 
A = 3 are M, = 0.6741 (40%) and M, = 0.8051 (60%). 

From Fig. [T^ we can see that for A = the perturbed 
configuration with its mass increased by 60% collapses 



-^ (A = 0), (//m) = 4.8 X 10"^ (A = 1) 
"^ (A = 2), (f/m) = 5.05 x 10"^ (A = 



FIG. 8. We show the Fourier transform of the evolution 
of the maximum value of the metric coefficient grr for the 
configuration with <j)-i — 0.28. The corresponding values of 
the quasi-normal frequencies for the different values of A are: 
(f/m) = 4.2 X 10" 
If/m) = 5.0 X 10~ 
3), while their fundamental frequencies are: (uj/m) = 0.91 
(A = 0), (tj/m) = 0.90 (A = 1), (a;/m) = 0.89 (A = 2) and 
[ui/m) = 0.88 (A = 3). Quasi-normal frequencies increase 
with A, a property that can also be noticed in Fig. [T] 



into a black hole, while for the rest of the A values this 
configuration migrates to another oscillaton located on 
the S-branch. The reason for this is that the A = 0- 
configuration has an initial mass that is larger than the 
critical mass of equilibrium configurations (see Table lll| . 
whereas for the other cases the initial mass is smaller than 
the critical one. The perturbed configurations with their 
mass increased by 40% are able to migrate to another 
S-oscillaton. 

For equilibrium configurations with (/'i(O) = 0.2, we 
increased its original mass in 25% obtaining the follow- 
ing initial masses: Mi = 0.6534 (A = 0), M, = 0.7198 
(A = 1), M, = 0.7919 (A = 2), and Mi = 0.8687 (A = 3). 
This configuration collapses into a black hole for all val- 
ues of A; in all cases, the initial mass of the perturbed 
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FIG. 9. The quasi-normal frequencies obtained from the evo- 
lution of slightly perturbed S-branch oscillatons as a function 
of the total mass Mt- Also shown is the migration of per- 
turbed S-branch and U-branch oscillatons, labelled A ~ E, 
see Figs. 



FIG. 10. We show the evolution of the total mass Mt and 
Rmax for different slightly-perturbed S-oscillatons up to a 
time t = 5000 m~^. S-oscillatons are not migrating nor de- 
caying, but only oscillating around their equilibrium positions 
in this plot. 



configuration is larger than critical one corresponding to 
each case. 

Finally, for equilibrium configurations with 0i(O) = 
0.44, we decreased its original mass by 20%: Mi = 0.4787 
(A = 0),M, = 0.5484 (A = 1), Mi = 0.6133 (A = 2), and 
Mi = 0.6696 (A = 3). As we can see in Fig. [H these 
configurations lose mass until they reach the position of 
another S-oscillaton. 

We have noticed that oscillatons maintain a fixed vi- 
bration frequency during its evolution. This is shown 
in Fig. I12[ where we present the migration path of os- 
cillatons with ^i(O) — 0.1 (with its mass increased by 
40%), and with </>i(0) = 0.44 (with its mass decreased 
by 20%), which are labeled A and B in Fig. [S] We can 
appreciate that the perturbed oscillaton with 0i(O) = 0.1 
(-1-40%) is migrating to an equilibrium configuration with 
01 (0) — 0.17. This can also can be seen in Fig. [131 which 
shows the profile of the metric coefficient §„ rapidly ap- 
proaching and oscillating around the final configuration. 

Then, as reported in Ref. 0], we have found that 
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FIG. 11. The scalar field <& and the radial metric function grr 
profiles of a strongly perturbed S-oscillatons with 0i(O) =0.1 
for the case with a mass increase of 40 percent. 
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FIG. 12. Evolution of different strongly perturbed S- 
oscillatons. If the initial mass is less than the critical value 
M^c., the configurations are able to evolve towards another 
S-configuration. Otherwise, the final fate may be the collapse 
into a black hole. 



FIG. 13. Evolved profiles of the metric coefficient Qrr for a 
strongly perturbed (/f>i(0) = 0.1-oscillaton for different values 
of A. Each plot indicates the final S-oscillaton each initial 
configuration is migrating to. 



strongly perturbed S-oscillaton are able to migrate to 
another S-oscillaton when their original mass is smaller 
than critical mass. But, if the original mass increases 
enough to be larger than the critical one, this perturbed 
configuration will collapse into a black hole, except in 
the case of diluted oscillatons, with low values of (^o, in 
which the collapse to a black hole can be prevent by the 
gravitational cooling mechanism. 



4- U -branch 

We identify as U-oscillatons the equilibrium configura- 
tions that are located on the right-hand side of the criti- 
cal configuration in a plot of Mt versus 4>i (0), see Fig. [21 
or located on the left-hand side in a plot of Mt versus 
Rmax, see Fig. [3l To evolve these equilibrium configu- 
rations, we also use the slightly-perturbed, by numerical 
inaccuracies, configurations, which in general decay and 
migrate to a configuration located on the S-branch. 



Fig. [14] shows some slightly perturbed U-oscillatons, 
and we can see that the equilibrium configuration is more 
unstable (the quicker it starts to migrate) for larger val- 
ues of (/'i(O). This also can be seen in Fig. [151 where 
we show the evolution of the maximum value of the 
radial metric coefficient g^ for the slightly perturbed 
configurations with central field values 4>i{0) — 0.5 and 
0i(O)=O.7. 

Thus, we confirm that the U-oscillatons are intrinsi- 
cally unstable under small perturbation, they decay and 
migrate to the S-branch. In Fig.[9l we show the migration 
of the slightly perturbed U-oscillaton with (j>i{Q) — 0.5 
and 01 (0) = 0.7, labeled C and S for A = 0, MO) = 0.7 
labeled E for A = 1 and A = 2, and 0i(O) = 0.58 labeled 
£: for A = 3. 



5. Perturbed U -branch 

We study now the behavior of the U-branch equilib- 
rium configurations under strong perturbations. An ex- 
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FIG. 14. Slightly-perturbed U-brancli configurations. In gen- 
eral, all of them migrate (left to right on the plots) to an 
equilibrium configuration located on the S-branch. 



FIG. 15. Evolve profiles of the maximum value of the radial 
metric function grr for the slightly perturbed configurations: 
MO) = 0.5 and (j>i{0) =0.7 



ample of the typical behaviors is provided by the config- 
uration with (pi = (0) = 0.6 for all values of A. First, 
we increase its original mass by 2%; the resulting initial 
masses are: Mi = 0.598 (A = 0), M^ = 0.678 (A = 1), 
M, = 0.734 (A = 2), and M^ = 0.755 (A = 3). 

In the case of A = 0, this configuration collapses into a 
black hole, as reported also in[3l- But, for other A values, 
the same configuration is able to migrate to the S-branch, 
see Fig. 1161 The migration path for this configuration is 
labeled as C, for A = 1,2,3, in Fig. [HI In contrast, a 
mass increase by 5% for the equilibrium configuration 
with 01 (0) = 0.5 provokes a rapid collapse into a black 
hole, independently of the value of A. 

In another experiment, we decrease the original mass 
of a 01 (0) = 0.6-oscillaton by 2%. The initial masses 
obtained are: Mi = 0.574 (A = 0), M^ = 0.651 (A == 1), 
M, = 0.706 (A = 2), and M, = 0.725 (A = 3). As 
expected, these perturbed U-oscillatons lose mass and 
migrate to the S-branch. The evolution of these strongly 
perturbed configurations appears in Fig. I16[ and their 
path migration is labeled as D in Fig. [S] for all the values 
of A. 



Fig. [17] shows the behavior of the metric functions Err- 
and gtt for the case of strongly perturbed U-oscillaton. 
In the plot for A = 0, we can see that gu shows the 
well known "collapse of the lapse", and g^ shows the 
" grid stretching" ; both phenomena signal the possible 
formation of a black hole. These behaviors are absent in 
the evolution of the metric functions for other values of 
A. 

Migration of the strongly-perturbed U-oscillaton cor- 
responding to 01 (0) = 0.6, with a mass decrease of 
2%, towards a S-oscillaton. See also Fig. [181 where we 
show the evolved profile of the radial metric function 
grr of the 01 (0) — 0.6-oscillaton migrating to the fol- 
lowing S-branch configurations: 0i(O) = 0.3 for A = 0, 
0i(O) = 0.28 for A = 1, 0i(O) = 0.18 for A = 2, and 
01 (0) =0.12 for A = 3. 



6. The S- U transition point 

We look carefully at the equilibrium configurations lo- 
cated nearby the critical configuration, .i.e., the most 
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FIG. 16. Evolution of the perturbed configurations witfi 
(/!>i(0) = 0.5 and 01(0) = 0.6 in a plot of Mt versus Rmax- 
Only in the case of A = 0, both configurations collapse into 
a black hole if their mass is increased, even by a little. For 
other values of A, the larger values of the critical mass are 
able to prevent the collapse of the less perturbed of the con- 
figurations. 



massive equilibrium configuration, see for instance Fig. [51 
which is usually called the S-U transition point 3] . Fig.fTOl 
shows the evolution of some slightly perturbed configu- 
ration near to the S-U transition point. With these ex- 
amples we confirm previous results: S-oscillatons are in- 
trinsically stable, whereas U-oscillatons are intrinsically 
unstable, no matter its proximity to the S-U transition 
point. Thus, this also confirms that the critical configu- 
ration is a true stability-instability transition point. 



V. CONCLUSIONS 

In this work, we solved numerically the Einstein-Klein- 
Gordon system for a massive and real scalar field en- 
dowed with a scalar potential containing a quartic self- 
interaction term. The diverse numerical experiments 
confirm the same general properties found in the case 



FIG. 17. Radial profiles of the metric functions grr and 
gtt for a strongly-perturbed U-configuration corresponding to 
(t>i{0) — 0.6 with an excess mass of 2%. The formation of a 
black hole seems to happen for A = 0, but its formation is 
prevented for other values. 



of the free massive case, which are also present in the 
case of boson stars. 

For all cases of the quartic self-interaction, it is found 
that there is a critical equilibrium configuration, which is 
the most massive one in each case; this critical configu- 
ration allows the separation of configurations into stable 
and unstable ones, the so-called S and U branches. As 
expected, the mass of the critical configuration increases 
for larger values of the quartic interaction. 

Equilibrium configurations located on the S-branch 
are intrinsically stable, and vibrate with definite quasi- 
normal frequencies when slightly perturbed. If strongly 
perturbed, they can either migrate to another S- 
configuration, or collapse into a black hole, if its initial 
mass is smaller or larger, respectively, than the critical 
one. On the other hand, U-oscillatons are intrinsically 
unstable configurations. Under small perturbations, they 
are able to lose mass and migrate towards and equilib- 
rium configuration on the S-branch, as long as its initial 
mass is not bigger than the critical one. However, un- 
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FIG. 18. The evolved profiles of the radial metric function 
Qrr of a perturbed U-oscillaton with (/!>i(0) — 0.6 and different 
values of A. In both cases, the evolution goes towards an 
equilibrium configuration in the S-branch. 



we speak of real or complex scalar fields. 
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FIG. 19. Evolution of the slightly perturbed configuration 
near the S-U transition point. The nature of the transition 
point does not depend upon the value of A. 



der strong perturbations, U-oscillatons cannot prevent 
its collapse into black holes. 

All in all, we have found that oscillatons with a quartic 
self-interaction share similar properties with their boson 
star counterparts. The larger the quartic interaction, the 
larger values the critical mass takes, but also the space 
for stable configurations is reduced because more diluted 
values of the scalar field are needed to sustain them. This 
can be seen in the shift of the critical scalar field value 
towards zero as the quartic interaction is increased. 

It proved quite difhcult to obtain equilibrium configu- 
rations for large values of the quartic interaction, because 
of the larger system of equations that must be solved in 
the case oscillatons (in terms of a Fourier expansion) as 
compared to the case of boson stars. Nonetheless, we 
have been able to explore the space parameter beyond the 
free case (A = 0), and provided evidence that points out 
to the close similarity of scalar configurations, whether 
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